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A CONVERGENT EXPLICIT FINITE DIFFERENCE SCHEME 
FOR A MECHANICAL MODEL FOR TUMOR GROWTH 

KONSTANTINA TRIVISA AND FRANZISKA WEBER 


Abstract. Mechanical models for tumor growth have been used extensively 
in recent years for the analysis of medical observations and for the prediction of 
cancer evolution based on imaging analysis. This work deals with the numer¬ 
ical approximation of a mechanical model for tumor growth and the analysis 
of its dynamics. The system under investigation is given by a multi-phase 
flow model; The densities of the different cells are governed by a transport 
equation for the evolution of tumor cells, whereas the velocity field is given 
by a Brinkman regularization of the classical Darcy’s law. An efficient finite 
difference scheme is proposed and shown to converge to a weak solution of the 
system. Our approach relies on convergence and compactness arguments in 
the spirit of Lions [23] . 


1. Introduction 

1.1. Motivation. Mechanical models for tumor growth are used extensively in re¬ 
cent years for the prediction of cancer evolution based on imaging analysis. Such 
models are based on the assumption that the growth of the tumor is mainly lim¬ 
ited by the competition for space. Mathematical modeling, analysis and numerical 
simulations together with experimental and clinical observations are essential com¬ 
ponents in the effort to enhance our understanding of the cancer development. The 
goal of this article is to make a further step in the investigation of such models by 
presenting a convergent explicit finite difference scheme for the numerical approxi¬ 
mation of a Hele-Shaw-type model for tumor growth and by providing its detailed 
mathematical analysis. Even though the main focus in the present work is on the 
investigation of the evolution of the proliferating cells, it provides a mathematical 
framework that can potentially accommodate more complex systems that account 
for the presence of nutrient and drug application. This will be the subject of future 
investigation [30]. 

1.2. Governing equations. In the present context the tissue is considered as a 
multi-phase fluid and the ability of the tumor to expand into a host tissue is then 
primarily driven by the cell division rate which depends on the local cell density 
and the mechanical pressure in the tumor. 

1.2.1. Transport equations for the evolution of the cell densities. The dynamics of 
the cell population density n{t, x) under pressure forces and cell multiplication is 
described by a transport equation 

dtn — diY{nu) = nG{p), x G t>0 (1.1) 

where n represents the number density of tumor cells, u the velocity field and p the 
pressure of the tumor, fl is a bounded domain in d = 2,3. The pressure law is 
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given by 

p{n) = ain?, (1-2) 

where 7 > 2. Following [3, 29], we assume that growth is directly related to the 
pressure through a function G(-) which satisfies 

GgC'I(K), G'(-)<-/?< 0 , G{Pm) = 0 for some Pm > 0 . (1.3) 

The pressure Pm is usually called homeostatic pressure. Here, and in what follows, 
for simplicity we let 

G{p)=a-Pp\ (1.4) 

for some a, 13,9 > 0 . 

1.2.2. The tumor tissue as a porous medium. The continuous motion of cells within 
the tumor region, typically due to proliferation, is represented by the velocity field 
u := VVF given by an alternative to Darcy’s equation known as Brinkman’s equation 

p = W- pAW (1.5) 

where /r is a positive constant describing the viscous like properties of tumor cells 
and p is the pressure given by ( 1 . 2 ). 

Relation (1.5) consists of two terms. The first term is the usual Darcy’s law, 
which in the present setting describes the tendency of cells to move down pressure 
gradients and results from the friction of the tumor cells with the extracellular ma¬ 
trix. The second term, on the other hand, is a dissipative force density (analogous 
to the Laplacian term that appears in the Navier-Stokes equation) and results from 
the internal cell friction due to cell volume changes. A second interpretation of 
relation (1.5) is the tumor tissue can be viewed as “fluid like.” In other words, the 
tumor cells flow through the fixed extracellular matrix like a flow through a porous 
medium, obeying Brinkman’s law. 

The resulting model, governed by the transport equation (1.1) for the population 
density of cells, the elliptic equation (1.5) for the velocity field and a state equation 
for the pressure law ( 1 . 2 ), now reads 

f dtn — div(nVVF) = an — , x G it, t > 0 

\ -pAW + W = an^. 

We complete the system (1.6) with a family of initial data no satisfying (for some 
constant C) 

no > 0 , p{no) < Pm, ||no||Li(Rd) < C. (1.7) 

The objective of this work is to establish the global existence of weak solutions 
to the nonlinear model for tumor growth ( 1 . 6 ) by designing an efficient numerical 
scheme for its approximation and by showing that this scheme converges when the 
mesh is refined. The main ingredients of our approach and contribution to the 
existing theory include: 

• The introduction of a suitable notion of solutions to the nonlinear system (1.6) 
consisting of the transport equation (1.1) and the Brinkman regularization 

(1-5). 

• The construction of an approximating procedure which relies on an artificial 
vanishing viscosity approximation and the establishment of the suitable com¬ 
pactness in order to pass into the limit and to conclude convergence to the 
original system (cf. Section 3, Lemma 3.7). 

• The design of an efficient numerical scheme for the numerical approximation 
of the nonlinear system (1.1)-(1.5). 


ON A TUMOR GROWTH MODEL 


3 


• The proof of the convergence of the numerical scheme. In the center of the 
analysis lies the proof of the strong convergence of the cell densities. This is 
achieved by establishing the weak continuity of the effective viscous pressure 
in the spirit of Lions [23] (cf. Section 4, Lemma 4.8). 

• The design of numerical experiments in order to establish that the finite differ¬ 
ence scheme is effective in computing approximate solutions to the nonlinear 
system (1.6) (cf. Section 5). 

For relevant results on the analysis and the numerical approximation of a two- 
phase flow model in porous media we refer the reader to [ 6 ]. Related results on the 
numerical approximation of compressible fluids employing the weak compactness 
tools developed by of Lions [23] in the discrete setting have been established by 
Karper et al. [19, 16, 17, 18] and Gallouet et al. [13]. 

Relevant work on the mathematical analysis of mechanical models of Hele-Shaw- 
type have been presented by Perthame et al. [25, 26, 27, 28]. The analysis in 
[27] establishes the existence of traveling wave solutions of the Hele-Shaw model 
of tumor growth with nutrient and presents numerical observations in two space 
dimensions. The present article is according to our knowledge the first article pre¬ 
senting rigorous analytical results on the global existence of general weak solutions 
to Hele-Shaw-type systems. 

A different approach yielding results on the global existence of weak solutions to 
a nonlinear model for tumor growth in a general moving domain fit C without 
any symmetry assumption and for finite large initial data is presented in [10, 8 , 9]. 
But in contrast to the present nonlinear system, the transport equation for the 
evolution of cancerous cells in [10, 9] has a source term which is linear with respect 
to cell density. 

Relevant results on nonlinear models for tumor growth governed by the Darcy’s 
law for the evolution of the velocity field are presented by Zhao [31] based on the 
farmework introduced by Friedman et al. [14, 5]. 

1.3. Outline. The paper is organized as follows: Section 1 presents the motivation, 
modeling and introduces the necessary preliminary material. Section 2 provides a 
weak formulation of the problem and states the main result. Section 3 is devoted 
to the global existence of solutions via a vanishing viscosity approximation. In 
Section 4 we present an efficient finite difference scheme for the approximation of 
the weak solution to system (1.6) on rectangular domains and Section 5 is devoted 
to numerical experiments. A discretized Aubin-Lions lemma and some technical 
lemmas are presented in Appendices A and B respectively. 


2. Weak formulation and main results 

Notation 2.1. For cp : (0,T) x D —>■ K, tp : (0,T) x D —>■ M.‘^, we will denote by 
Vyj := S/xP = ..., and div ip := diva, cp = Y.i=i the gradient 

and divergence in the spatial direction in fl. 

2.1. Weak solutions. 

Definition 2.2. Let fl a bounded domain in d = 2,3, which is either rectan¬ 
gular or has a smooth boundary dfl and T > 0 a finite time horizon. We say that 
{n,W,p) is a weak solution of problem (1.1)-(1.5) supplemented with initial data 
{no,Wo,po) satisfying (1.7) provided that the following hold: 

• {n,W,p) > 0 represents a weak solution of (1.1)-(1.5) on (0,r) x O, i.e., for 
any test function p S C(?°([0,r] x IR‘^),T > 0, the following integral relations hold 
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/B'i 


nip^T, ■) dx 


In particular, 


- no(p{0,-)dx 

Jm.‘‘ 



n\/W ■ V(/3 + nG{p)(p{t, •)) dxdt. 

( 2 . 1 ) 


ne LP((0,T) X n), for all p> 1. 


We remark that in the weak formulation, it is convenient that the equations (1.1) 
hold in the whole space provided that the densities n are extended to be zero 
outside the tumor domain. 


• Brinkman’s equation (1.5) holds in the sense of distributions, i.e., for any test 
function p G C'“(M‘’*) satisfying 


(f\oa = 0 for any t G [0,r], 
the following integral relation holds for a.e. t G [0,r], 


J arpLpdx = J ^/iVW . V(/3 + 


dx. 


( 2 . 2 ) 


and p = almost everywhere. All quantities in (2.2) are required to be integrable, 
and in particular, W G L°°([0,T]; 


The main result of the article now follows. 


Theorem 2.3. Let fi C be a bounded domain with smooth boundary dLl, 0 < 
T < oo. Assume that the initial data uq G L°°(Ll) with 0 < uq < n^o '■= Pm'^ o,nd 
that G{-) is of the form (1.4). Then the problem (1.1)-(1.5), admits a weak solution 
in the sense specified in Definition 2.2. 

The following two remarks are now in order. 

Remark 2.4. In Section 3, such a solution is obtained as the limit of the vanishing 
viscosity approximations {ne,We,Pe) of (3.1) to (1.6) as £ —)■ 0. 

Remark 2.5. In Section 4, such a solution is obtained in the case of a rectangular 
domain, as the limit of the sequence of approximations {nh,Wh,Ph) computed by 
the numerical scheme (4.1) - (4.3) as h —>■ 0. 


3. Global existence via vanishing viscosity 

In this section we prove Theorem 2.3 by constructing an approximating scheme 
which relies on the addition of an artihcial vanishing viscosity approximation 

{ dtUe — div(n£VITe) = an^ — + EAn^, x € ft, t >0 

pAWe -We = an], (3.1) 

ne( 0 ,.) = ng, 

where ng is a smoothnened version of no, that is ng = no * Te lor a smooth function 
(fg with compact support, and a bounded domain kl G K'’* with smooth boundary 
or alternatively the d-dimensional torus T'^, and we establish its convergence to the 
nonlinear system (1.6) at the continuous level. For simplicity, we assume a = 1 
and homogeneous Neumann boundary conditions for n^ and Wg (if the domain is 
a torus we can also use periodic boundary conditions). 

Theorem 3.1. For every £ > 0, the parabolic-elliptic system (3.1) admits a unique 
smooth solution (jig,Wg,Pg). 
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Proof. The proof of this result relies on classical arguments (cf. Ladyzhenskaya 
[20]), namely by employing the Contraction Mapping Principle and the regularity 
of the initial data one can show the existence of a unique solution ing,Wg,pg) 
defined for a small time T > 0. Then one derives apriori estimates establishing 
that the solution does not blow up and in fact is defined for every time. Finally, a 
bootstrap argument yields the smoothness of the solution. □ 

The remaining part of this section aims to establish the necessary compactness 
of the approximate sequence of solutions {rig, We,pg). 


3.1. A priori estimates. We start by proving that rig are uniformly bounded 
independent of e > 0 and nonnegative: 

Lemma 3.2. If 0 < ng{0,-) < rioo := Pm'* < oo uniformly in e > Q, then for 
any < > 0, the functions ng{t, •) are uniformly {in e > 0) bounded and nonnegative, 
specifically, 

0 < Tahing{t, x) < ma.xng{t, x) < rioo- 

{t,x) 

Proof. First we notice that if Wg has a maximum at a point xq, then AWe(-, xq ) < 0 
and therefore Wg = pg + pAWg < pg. Similarly, if it has a minimum at a point xq, it 
will satisfy AWg{-,XQ) > 0 and therefore Wg > pg. If Wg attains a strict maximum 
on the boundary, i.e., there is a point Xq S dfl such that Wg{xo) > Wg{x) for 
any other cc S fl, we apply Hopf’s Lemma, e.g. [12, p. 347], to the function 
V := Wg — maxp 2:) Pe{t, x) which satisfies 

—fiAv + V = Pg — maxpg{t, x) < 0, 

(t,x) 

which has a strict maximum at the point xq. If v{xa) < 0, then Wg < Wg{xo) < 
maxp_2,) a;) and otherwise Hopf lemma gives VWg{xo) ■ v = Vv{xo) • u > 0 
where we have denoted the boundary normal u, this contradicts the homogeneous 
boundary conditions. In a similar way we show that Wg > min(t a,) pg{t, x) (applying 
Hopf’s lemma to —Wg and hence 


m.m.pg{t,x) < Wg < maxpg{t,x). 

{t,x) 


(3.2) 


We rewrite the evolution equation for Ug using the equation for the potential Wg, 

dtUg - VWg ■ Vug = ngG{pg) + -ng{pg - Wg) + eAug. (3.3) 

Now assume {to,xo) is a point, where ng{to,Xo) > Uoo reaches its maximum (and 
therefore also pg{to,xo) > Pm reaches a maximum). Then 'Vng{to,Xo) = 0 and 
Ang{to,xo) < 0. Hence 

dtng{to,Xo) < ngG{pg) + -ng{pg - Wg). 

By (3.2), the second term on the right hand side is nonpositive and since G{pg{tQ, Xq)) < 
0 for Pg > Pm, we get 

dtng{to,xo) < 0 . 

Hence Ug will decrease and if initially Rq A Roo this implies that ng{t, •) < Rco for 
any later time t > 0. To show the nonnegativity of Ug, we integrate the evolution 
equation for Ug, 


^ [ Ugdx = [ ngG{pg)dx. 

“t Jn Jo. 
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On the other hand, multiplying the same equation by a regularized version of the 
sign function, integrating and then passing to the limit in the approximation, we 
have 

d 


dt. 


\ne\dx < / \ni;\G{ps)dx, 


Subtracting the two equations from one another, and using that \ne\ — > 0, 

4 / \\ne\ - ne\dx < [ \\ns\- nJG{pe)dx, 

Jq Jq. 

< max |G(s)| / \\ne\ — rLAdx. 
sG[0,Pm] jo 

Now using Gronwall’s inequality and that |no| ~ pq = 0 by assumption, we obtain 

/ I IueI — Tie I {t)dx = 0 

JO 

and thus that x) > 0 almost everywhere. □ 

Next we prove a simple lemma on the regularity of Wg. 

Lemma 3.3. We have that 

Wg C 

for any q G [1, oo) uniformly in e > 0 and 

Wg,mig C L~((0,T) X O)), 

uniformly in e > 0 as well. 

Proof. We square the equation for Wg and integrate it over the spatial domain and 
then use integration by parts, 

[ \pg\'^dx= [ \Wg\^-2p,WgAWg + p^\AWg\^dx 
Jn Jq 

= f IWgl"^ + 2p\VWg\^ + p^\V^Wg\‘^dx. 

Jq 

By the previous Lemma 3.2, we have that pg is uniformly bounded in e > 0 and 
therefore that the left hand side of the above equation is bounded and that Wg € 
L°°([0,Tj; Using a Calderon-Zygmund inequality (e.g. [15, Thm. 9.11.]), 

we obtain Wg G L°°{[Q,T]]W'^'^{Wj) for all q G [l,oo). By the Sobolev embedding 
theorem, this implies that in particular VWg G L°°{{0,T) x fl). The second claim 
follows from (3.2) and the uniform bound on the pressure proved in Lemma 3.2. □ 

3.2. Entropy inequalities for Ug. To prove strong convergence of the approxi¬ 
mating sequence {{ng,Wg,pg)}g:^o, it will be useful to derive entropy inequalities 
for ng. To this end, the following lemma will be useful: 

Lemma 3.4. Let / : M —>■ K fee a smooth convex, nonnegative function and denote 
fg := fijig). Then fg satisfies the following identity 


dtfe - dW{fgVWg) - eAfing) 

= {f'{'ng)ng - fg)AWg + f{ng)ngG{pg) - ef"{ng)\Vng\'^ 

where 


f [ f"{ng)\Vngfdxdt<C, 
Jo Jq 


(3.4) 


(3.5) 


with C > 0 a constant independent of e > 0. In particular, this implies that 
dtfe = 9s + kg with gg G L^([0,r] x LI) and kg G L^{[0,T]; 1U“^’^(0)). 
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Proof. The identity (3.4) follows after multiplying the evolution equation for Ug, 
(3.3), by /'(re) and using chain rule. Integrating the inequality in space and time, 
we obtain 


fe{T)dx + e f f 

1 Jo Jn 


/"(RE)|VnEp dxdt 


= [ fA0)dx+ [ [ (/'(r 


fe)^Ws + f dxdt 


Jn Jo Jn 

The right hand side is bounded by the assumptions on the initial data and the L°°- 
bounds proved in Lemmas 3.2 and 3.3. This implies (3.5). Therefore the right hand 
side of (3.4) is contained in L^((0, T) x ^l). Using (3.5) for the third term on the left 
hand side, we conclude that it is contained in L^{[0,T]-, The second term 

on the left hand side is contained in L°° Hence dtfe = Qe + 
with G L^([ 0 ,T] x H) and fcg G L^([ 0 , T]; IU“^’^(fl)) and in particular, dtfe G 
(H)) by the Sobolev embedding (1* = d/{d— 1)). □ 


Remark 3.5. The preceeding lemma implies that the time derivative of the approx¬ 
imation of the pressure dtPe = dt\nhp = Ps + ke where pe is uniformly bounded 
in L^{[0,T] X fl) and in L^{[0,T]] Hence dtWe = Ug + Ve where 

Ue G L\[0,T];H^{n)) solves -pAUe + Ue = k^ and G L^([0,T];W^’^(fl)), 
1 < r < 1* solves —pAVe + Ve = Pe (see [1, Thm. 6.1] for a proof of the second 
statement). Hence dtWe G L^([0,r]; IU^’’’(r2)) for any 1 < r < 1*. 

3.3. Passing to the limit e —)■ 0. The estimates of the previous (sub)sections 
allow us to pass to the limit e —)■ 0 in a subsequence, still denoted e, and conclude 
the existence of limit functions 


tt-e ^ r > 0, inL'^([0,T] x H), 1 < g < oo, 
Pe^P>0, inL'^([0,T] X fl), 1 < g < 00, 


where Pe := n] and 0 < n,p G L°°([0,T] x H). Using Aubin-Lions’ lemma for Wg 
and VWe, we obtain strong convergence of a subsequence in L'^{[0,T] x H) for any 
g G [0, oo) to limit functions W, VIU G L'^{[0, T] x H). Moreover, from the estimates 
in Lemma 3.3 we obtain that W G L°°([0, T] x fl) n L°°([0, T]; 1U^’'^(H)). Hence we 
have that {n,W,p) satisfy for any € C'q([0,T) x H), 


n mpt — n'VW ■ \/(p dxdt + / tiq plfi, x)dx = — / / nG{p)(pdxdt 

! Jn Jo Jn 


Wtjj + AtVIU • \/ip dxdt = 


0 Jn 


(3.6) 


p tf dxdt 


0 Jn 


where nG{p) is the weak limit of neG{pe). To conclude that the limit {n,W,p) is 
a weak solution of ( 1 . 6 ), we need to show that rig converges strongly and therefore 
in the limit p = p -.= rp and nG{p) = nG{p). For this purpose, we combine a 
compensated compactness property (Lemma 3.7) with a monotonicity argument. 
We will also make use of the following lemma which was proved in a more general 
version in [7, 24]; 


Lemma 3.6. Let ri, f G L°°([0,T] x H) and u G L°°([0, T]; iJ^(H)) with divM G 
L°°([0,T] X H) satisfy 

nt - div(Mn) = /, (3.7) 

in the sense of distributions. Then for all continuously differentiable functions 

6g C'i(K), 

b{n)t — div{ub{n)) = b'{n)f + \b'{n)n — b{n)] divii, 
in the sense of distributions. 


(3.8) 
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Proof. We let 0 < ^ G be a smooth, radially symmetric mollifier, i.e. 

ipix) = %f{—x) and J^a+i 'f’{x)dx, with supp('(/)) C -Bi(O) and denote for <5 > 0, 
tps{x) '■= S~^'^~^^'>^|:{x/6). Then we choose as a test function in (3.7) 'ijjs{s,y)ip{t + 
s,x + y), with If is compactly supported in {6,T — S) x where includes all 
the points a; in which have distance d{x, dfl) > S and do a change of variables: 


n n{t -s,x- y)'ips{s, y)dtip{t, x) -n{t-s,x- y)u{t, x)^fs{s, y) ■ 'V(p{t, x) dxdt 
! 

n f{t - s,x- ?/)' 05 (s, y)tf{t, x) dxdt. 

! 

Integrating in {s,y), this becomes 

n (n * 'il)i){t,x)dty}{t,x) — {nu) * ifsit^x) ■ Vip(t,x) dxdt 
! 

= - / / if *ips){t,x)(f{t,x)dxdt. 

Jo Jn 

We define ns := n* tps and fs:=f* ips and choose as a test function ip := b'{ns)(j) 
for a smooth (p compactly supported in {6, T — S) x (which is possible since ns 
is smooth and bounded thanks to the convolution.). Then we can rewrite the last 
identity using chain rule as 

n b{ns)dt4> — b{ns)u ■ \7cj)dxdt 
i 

= - / {b'{ns)fs + [b'{ns)ns-b{ns)]divu + b'{ns)rs)(j)dxdt. 

Jo Jq 

where rs := div((mx) *ips) — div(n5M). By [22, Lemma 2.3], we have that —>■ 0 in 

X fl) and thanks to the properties of the convolution that b{ns) —>■ b{n) 
almost everywhere as well as fs ^ f a.e. when i5 —> 0. Thus we obtain that in the 
limit <5 —> 0, n satisfies 

f f b{n)dt4i — b{n)u ■ \7(j)dxdt = — f f {b'{n)f + [b'{n)n — b{n)] div u) p dxdt. 

Jo Jn Jo Jn 

which is exactly (3.8) in the sense of distributions. □ 

Applying Lemma 3.6 for the weak limit n in (3.6) with b{n) = n^, we obtain 
that n satisfies 

n n^ipt — -Vpdxdt = — f f {2nnG{p) + n"^ AW)p dxdt (3.9) 

! Jo Jn 

for any test functions p G C'q((0,T) x fl). On the other hand, from (3.4) for 
b(n) = we obtain after integrating in space and time 

f nl(T)dx— f nl{0)dx< f f n^AW^ + 2nlG{pe) dxdt 

Jq Jq, Jo Jq, 


Passing to the limit e —> 0 in this inequality, we have 

f n^{T) dx— 


In 


In 


riQ dx < 


lo Jn 


n'^AW + 2n^G{p) dxdt, 


(3.10) 


where denotes the weak limit of and n'^AW and n‘^G{p) are the weak limits 
of n^AWe and nlG{pe) respectively. Letting r —>■ 0 in this inequality, we obtain, 
thanks to the boundedness of the integrand on the right hand side. 


/ n^(0)fi2:— / Rq da; < 0. 

Jn Jn 
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On the other hand, since b{n) = is convex, we have and hence n^(0, x) = 


nl{x). 


We now choose smooth test functions (p^ approximating ip(t, x) = 1[o,t] (t), where 
T G (0,T'], in inequality (3.9) and then pass to the limit in the approximation to 
obtain the inequality 

f v?‘{T)dx— f nldx= f ( {2nnG{p) + n^AW) dxdt (3-11) 
t/Q '/Q J 0 J ^ 

Subtracting (3.11) from (3.10), we have 

J (n^ — {T)dx 

<JJ (2n'^G(p) — 2nnG{p) + AW + n'^AW — n'^AW^ dx dt. 

(3.12) 

Now using the explicit expression of G, (1.4), the first term on the right hand side 
can be estimated as follows: 

J J (^2n^G{p) — 2nnG{p)j dx dt 

= 2 / fa — P ^71^+7® — 7171^+^®^ dxdt 

0 r2 

/ / a ^77^ — 77^^ — P — 77^+9'®^ dx dt 

»/ 0 


(3.13) 


< 2 

< 2q; 


77^ — 77^ ) dx dt 


) 


where we have used [24, Lemma 3.35], which implies , for the first 

inequality. To estimate the second term on the right hand side, we use that AW 
is bounded thanks to Lemma 3.3 and that > nP by the convexity of f{x) = x'^. 
Hence ^ ^ 

[ [ AW (vP — vP\ dxdt < —^ [ [ (vP — 77^') dxdt. (3.14) 

Jo Jn ^ Jo ' 

For the last term, we use the following lemma. 

Lemma 3.7. The weak limits {n,W,p) of the sequences {(tt^, lFe,Pe)}E>o satisfy 
for smooth functions S' : K —>■ K, 

J (s(77 )A1F - S^Awj dx = i y (p SH - pS{n)^ dx (3.15) 


where S{n)AW, S{n), pS(n) are the weak limits of S{ne)AWs, S{ne) and PeSijie) 
respectively. 

Applying this lemma to the second term in (3.12) with S{n) = 77 ^, we can 
estimate it by 

1 ( 


10 Jn 


77^A1F — 77^ AIF ) dx = — J yp n'^ — pn^j dxdt 

^ ^ 77 T' 77^ — 77^+^ ) dxdt 


M Jn 


J 


< 0 , 


using that (cf. [24]). Thus, 
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Hence Gronwall’s inequality implies 

IT? — (r)(ix < 0 

By convexity of the function f{x) = we also have < rfi almost everywhere 
and so 

n^{t, x) = x) 

almost everywhere in (0, T) x H. Therefore we conclude that the functions rig 
converge strongly to n almost everywhere and in particular also p = ri^ which 
means that the limit (n, W,p) is a weak solution of the equations (1.6). 



Proof of Lemma 3.7. We multiply the equation for Wg by S{ng) and integrate over 

n, 

[ fj,AWgS{ng) - WgS{ng) dx = - [ pgS{ng) dx. 

Jn Jn 

Passing to the limit e —>■ 0, we obtain 

f pAWS{n) — WS{n) dx = — ( pS{n)dx. (3.16) 

Jq Jq 

On the other hand, using the smooth function S{ng) as a test function in the weak 
formulation of the limit equation 

-pAW + W =p, 

and passing to the limit e —> 0, we obtain 


/ fiAWS{n) — WS{n) dx = — / pS{n)dx. 
JQ JQ 

Combining the last identity with (3.16), we obtain (3.15). 


□ 


4. Global existence via a numerical approximation 


We consider the problem in two space dimensions in a rectangular domain, for 
simplicity we use fl = [0,1]^, the generalization to other rectangular domains as 
well as three space dimensions is straightforward but more cumbersome in terms 
of notation, for this reason we restrict ourself to a square two dimensional domain 
here. For simplicity, we will also assume a = 1 in the Brinkman law in (1.6). We let 
h > 0 the mesh width, and At the time step size. We will determine the necessary 
ratio between h and At later on. For i,j = 1, ..., N^, where = 1/h, h chosen 
such that Nx is an integer, we denote grid cells Cij := {{i — l)h,ih] x {{j — l)h,jh] 
with cell midpoints Xij = {{i — l/2)h, (j —l/2)/i). In addition, we denote P” = mAt, 
m = 0,... Nt, where Nt = T/At for some final time T > 0. The approximation of 
a function / at grid point Xij and time t™ will be denoted ffj. We also introduce 
the finite differences, 


X'l Jij — 


Dff^j = ± 


fij 


Dff^ = ± 


j^m±L _ jtr, 


and define the discrete Laplacian, divergence and gradient operators based on these, 


:= {D^,DfY, div± + Dtf\% A, := div± 


For ease of notation, we also let and Vijjgi /2 denote the discrete velocities 

in the transport equation, specifically, given Wij, we let 


(4.1) 
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4.1. An explicit finite difference scheme. Given at time step m, 

we define the quantities at the next time step by 




pZ ■■= Ki-r- 


D+nZ + = nZG{pZ), 

where pij = (nij )’’' and the fluxes j = 1,2 are defined by 

„m I m 7 

l(l) / m _ m by b+i,j n-, , „+ „ 


pO; m'l _ _,,m 

'^i+l/ 2 yv^ )— ^i+l/ 2 ,j 


- i^K+i/2ZDI'nZj 


F 


( 2 ) 




) = 


iy+i/ 2 V‘^ O" ) --'Ji,j+ii 2 - 2 

We use homogeneous Neumann or periodic boundary conditions for both variables: 


b.i+i 


- Zvi,j+i/2\Dtn^j. 


(4.2a) 

(4.2b) 

(4.2c) 


(4.3) 




^rn _ m 

J = I 7 ■ 


<0 

— ^i,li 

m _ 171 

i = 1 ,. 



= IT,’”- 

TJ/m _ jj/m 

J = I 7 ■ 

• ■ 7 Nj; , 


= 

wm _ 

i = 1 ,. 

..,A,. 


The initial condition we approximate taking averages over the cells, 

Kj = i7^ / ^^ 0 ( 2 ;) dx, plj = \nljp, i,j = 1,..., No,. 

No I JCij 


4.2. Estimates on approximations. In the following, we will prove estimates 
on the discrete quantities obtained using the scheme (4.1)-(4.3). We 

therefore define the piecewise constant functions 

Nt 

fh{t,x) = ^ ^ fZ lc.,(a:)l[t-,t-+i)(t), {t,x) G [0,T] x ff, (4.4) 

m —0 iyj —1 

where / G {n,W,p'\. We first prove that stays nonnegative and uniformly 
bounded from above. 


Lemma 4.1. If 0 < < riao ■= < 00 uniformly in h > 0 and the timestep 

At satisfies the CFL condition 


At < min 


_ h _ p \ 

Smaxy \VhWff^\ + hG-° ’ J 


(4.5) 


(where G°° := maxsg]a+ G(s)), then for any t > 0, the functions nh(t, •) are uni¬ 
formly (in h > 0) bounded and nonnegative, specifically, defining rioo = Roo + 
4Atsup,>o (sN7G(s)), we have for all m> 0, 

0 < min n™ < max n™,- < ffoo • 

i,j i,j 


Proof. The proof goes by induction on the timestep m. Clearly, by the assumptions, 
we have 0 < < rioo- For the induction step we therefore assume that this holds 

for timestep m > 0 and show that it implies the nonnegativity and boundedness at 
timestep m + 1 . 

We first show that the Wfj are bounded in terms of the p™ . To do so, let us 
assume it has a local maximum Wft in a cell Cy, for some i, j G {!,..., N^}. Then 

DtWF<0, -D^WP1<0, k = l,2, 
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(if i or j G then because of the Neumann boundary conditions, the for¬ 

ward/backward difference in direction of the boundary is zero and thus the previous 
inequality is true as well). Hence 

= X E < 0. 

fc=i 


Therefore, 


1 


= Ci + < Pts < “ax In™.r. 

Similarly, at a local minimum WT) of Wh, we have 


and hence 


which implies 


Thus, 


> 0, -D-Wr, >0, fc = 1, 2 


> 0 , 


= Ci + > Ci > ““ K, r > 0- 


0 < W/t < max |n™ |'''. 

i,j 


(4.6) 


Now we rewrite the scheme (4.2c) as 


n™+^ = 
'bj 


+ af}^) nT,, +/3™n^i., +C,Z<-i,j +C.<.+i (4-7) 


where 


>ir=i-K[(K 


,) 


2/j LM“'*+l/2,jl ■*" “i+l/2,jl ll“i-l/2,jl “i-l/2,jJ 

+ (l<, + l/ 2 l + <, + 1/2) + (K,-l/ 2 l - <,-1/2)] 


(2).m 


= AtG(p™)+^ 


h . 


^2+1/2,j "T ^2,j+l/2 ^2 j'-1/2 


A” = ^ (<+ 1 / 2 ,, + i»r+./2,ji) 

c = |)(i«r-i/2,,i-"r-./2.i) 

^ At 


2h 

fl™ — 


(< 


i+ 1/2 




*j-l/2l ‘^i,i-l/2 


We note that P^jt Cij j Vi^j t ^ 0, and that under the CFL-condition (4.5), also 
“1^’™ + Hence, assuming that n™ > 0 for all i,j, we have 

m+1 \ / om I /-m , , nm \ \ 


( 


(l),m , (2),m 

a* j + 0^1' 


* j 


> 0 . 


We proceed to showing the boundedness of n^. Thanks to the CFL-condition (4.5), 
we have 
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Moreover, alj = 1- Using the induction hypothesis that 
n™- < fioo for all i,j and the nonnegativity of nn which we have just proved, we 
can estimate 






+ a. 


(2),m 


I f om I /-m , ra , /nm — 


= Tlao — 2 (rico — ^ij ) U 


(4.8) 


(2),m m 

J' '■jj' 


We can rewrite and bound using the equation for W^j, (4.2a), 

= At (G(p™ ) + 

= At (^g(p-) + ^(m^;:)-c.)) 

< At (^G(p-) + ^ (^2o - K-r)) 


< At ( G(p-) + {n^ - n-) 




1 


< AtG(p™ ) + — (n^ - n-), 

where we have used (4.6) for the first inequality, that /(a) — /(6) = f'(a){a — b) for 
some intermediate value d G [6, a], with f(a) = , for the second inequality and 

the CFL-condition for the last inequality. Now going back to (4.8) and inserting 
this there, we obtain, 




<n^-l{n^- <,) + (AtG(p™ ) + ^(ft^ - n™ ) 


<5noo + Jn- +Atn-G(p™) 


(4.9) 


If n'^j > «oo then G{p'!^j) < 0 and hence the expression in (4.9) is bounded by rioo- 
On the other hand, if n™- < n-oo, we can bound it by 


m+l ^ ±- 
^^,0 - 4*^° 

3_ 

- 4 ”'°° 


-n™ + Atn™.G(p™ ) 

7 fnoo+4Atsup (s^Ag(s)') 
4 V .«>n V / 


= no 


where we used the definition of Uao for the last equality. This proves that < 

rioo for all i,j if the same holds already for the n™-. □ 

Remark 4.2. The estimates in the proof of the previous lemma are very coarse and 
therefore one can use a much larger CFL-condition than (4.5) in practice. Also 
note that riao —t riao when At —> 0. 

4.2.1. Estimates on the discrete potential Wh- 

Lemma 4.3. We have that 

Wh,VhWh,VlWh C L°°i[0,T];L^in)), 

uniformly in h > 0, where Vh ■= 

Wh,AhWhCL^i{0,T)xn)), 
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uniformly in h > 0 as well. 

Proof. To obtain the L^-estimates, we square the equation for the potential Wh, 
(4.2a) and sum over all i,j, 

- 2n £ T 

i,j = l i,j=i iyj^l i,j=i 

Using summation by parts and that W satisfies either periodic or homogeneous 
Neumann boundary conditions, we obtain 

p" E +2/^ E + E = E 

i,j—i ^ii=i 

From the previous estimates, we know that nt S L°°{[0,T] x 17) uniformly in 
h > 0 and therefore also uniformly bounded in any other L^-space, which implies 
together with the above identity, that Wh,VhWh,'^hWh € L^([0,T] x 17). That 
Wh is uniformly bounded follows from (4.6) and the uniform bound on nu which 
was proved in the previous Lemma 4.1. 

Using this and the uniform boundedness of the pressure, we conclude by (4.2a) 
that also A^lU/j is uniformly bounded. □ 


Remark 4.4. Using the discrete Gagliardo-Nirenberg-Sobolev inequality, [2, Thm. 3.4], 
we obtain that VhWh G T°°([0, T]; L'J(17)) for 1 < g < q* = 2d/{d— 2). 

4.3. Discrete entropy inequalities for Uh- To prove strong convergence of the 
approximating sequence {(jih, W/t)}ii> 0 ; h will be useful to derive entropy inequal¬ 
ities for Uh- To this end, the following lemma will be useful: 


Lemma 4.5. Let / : K —>■ K fee a smooth convex function and assume that At 
satisfies the CFL-condition 


At < 




(4.10) 


I 16maxij |V;ilF™ | ’ Smax^ |Vft,lU™ | J 

Denote := f(nW) and fh a piecewise constant interpolation of it as in (4.4). 
Then /") satisfies the following identity 


Dtf- =\Df (<+1/2, im+mw)) +1^2 (<,+1/2 +<+i)) (4.11) 




+ ^D+ /'(nr.)|u™ 


,tl“z-l/2,,<l 


-rDZ 




/'«)K + l/2l« 


(4.12) 


/'«)K-1/2<2-< 


h^ 

Wo- 

4 2 


+ |2 


/ (’^i+l/2,)'*^i+l/2j'|41l 


r(j 




l« 


(4.13) 

(4.14) 

(4.15) 


- J/"«-i/2,)l<-i/2,ll«,f - J/"«-i/2)l<-i/2ll« 


— |2 
j' 


(4.16) 


- Tr(n((ti/2,)Ki/2,ii«,r - T/"«-+i/2)K+i/2iKrir 


+ (/'«X - + f'{nTMZG{pT,,) 


(4.17) 

(4.18) 
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+ ^f'in. 




)\Dr< 


(4.19) 


+i/2.i 




w/iere e [min{n™ ,max{n™ ,n™± 1 / 2 ,± 1/2 e 

[min{n™ ,n™-±J,max{n™-,n™±J] andnff^^'^ e [min{ri™ ,n^+^j,max{n™ ,ri((*+i}] 
and where the term (4.18) is uniformly bounded and the terms (4.16) - (4.17) and 
(4.19) satisfy 

hd+l/\f 

m—0 i,j 

hd+lAf 

E E/"K.+i/2)K,+i/2ll^2+<,f < O (4.20) 

m—0 i,j 
hd/\i2 d^T 

^^nKj )\DtnZ\^ < C, 

m—0 i,j 

In particular, this implies that the piecewise constant interpolation f^ is of the 
form fh = Qh + kh where gn e L^([0,r] x 11) and kh e L°°{[0,T];W~^'‘‘{il)) for 
any 1 < q < 00 if d = 2 and for 1 < q < q* = 2d/(d — 2) if d > 2, uniformly in 

/i > 0. 


Proof. We first rewrite the scheme for n™ as 


n+T7™ — r>+n^ 4- 

dd-ij — 2^i+l/2,j.^l + 2^d-l/2,j 


DfnZ 


-71™ r)+77™ -I- ^71™ r>~n'^ 

2^ij + l/2-^2 ”-7j + 2^ij-l/2^2 ddi^j 




l'^i+l/2j‘ 


\Dtd 


:D-, 


'^ij + 1/2 


\Dt< 


+ n™A,W-+n-G(p™). 

Then, using the Taylor expansion, 

f{b) - f{a) = f{a){b - a) + /"(a) 
where d € [min{a, b}, max{a, &}], we can write 


(a-5)2 


Dtm = f'Kj)DtnZ + 


-ji jpm 


m |2 

J' 


Dfff:, = /'(n-)Z4f n- ± 2 /" («Ili/2.,)|4?f <,r 

Dtf^j = f'K,)DUZ ± \f'i^T,,±i,2)\DtnTf 
D^fKj) = f\K±i/2,,)D^<j 
D^f\<j) = f\nT,j±i/2)D2nZ, 


(4.21) 


where n™ i/ 2 j> ’^rj±i/ 2 > ^Si/ 2 .j and ± 1/2 are intermediate values. Hence, 

multiplying equation (4.21) by f’in^j), it becomes 


Diff:, = ^f"Kr^'')\Dt<f 

1 ' rr , arm h 


+ ^m |2 


+ ^<+1/2,Ptfr:, - T/"(n'li/2.,)«)Ti/2j44fn™ 


4 

+ + jnK-i/2>T-i/2jDfn^^\^ 
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+ +1/2^2+- Jr(n™+i/2X,+i/2i^2+<, 


m |2 


4^ 

h 


— |2 


+ 2<,-1/2^2-/™ + 4/"(-r.-l/2)<,-l/2l^2-<, 

+ [/'«,)k+i/2.,i^i+<, 1 - 5/"(^r-i/2,,)kr-i/2,,ii^r< 


4 

h „- 
4 


— |2 


/'(n™ )K, + l/2l^2+<. - T/"(^r2-l/2)K.-l/2ll^2-^^ 


+ 4^r [r«)i^^r-i/2j^r<, - 4/"(nr+i/2,,)kr+i/2,,ii^r< 


4- 

h 


— |2 


h „+ 
4 


4- 


+ ^m |2 

4 I 


+ ^m |2 


r(n-)l<,-i/2l^2"<,J - ^nKj+i/2)Kj+i/2\\Dtnl 

+ /'(n™. X, + f{nZ)nZG{p^^) 

\d-, (xi/2. (/i:;-+/;?i4) + (<,+ 1/2 (/^+/^;+i)) 

/'«-)IXl/2j^f<J +X2" [/'«,)l<, + l/2l^2+< 

/'«,)l«r-l/2j^r<,l + X2^ f/'«,)K,-l/2l^2-< 


+ 2^ 
h 

+ 4^^" 
h . 

4 1 


rx 


.\Dt< 




.i + l/2^‘'j.i + l/2l^2 "-ij 


IX-- 12 


-X 

4 2 

- j/"X-l/2.,)IXl/2jlXXP - ^/"X-l/2)K,-l/2ll^2-Xf 

- \f'{K+l/2,)\<+l/2,\\DtnZ? - JrX + l/2)K, + l/2ll^2+< 


|2 _ 

4^ 

£f (^ra \^m , 


m |2 


+ (/'XX, - /IS)AX™ + /'XX-GX)- 

which implies (4.11)-(4.19). In particular, for /(x) = this becomes 

A+/™ = At|A+Xf 

+ lot X+/.-i.,)) + (<,- 1/2 im+kj-i)) 


- -^d: 


Xl/2,,lXXf 


/i2 

- 


Xl/2l^: 


+ ^m |2 


XIXl/2.,lXX 

h 

+ -X 

2 ^ 

XIXl/2lX</ 

XIXi/2.,IXX' 

h . 

+ -X 

2 ^ 

x>r,-i/2ixx 


- Jx-i/2.,lixxf - ^IXl/2llXXf 

- JlX/2,,llXXf - JlXi/2llX< 

+ /”!-A,X + 2/”!-GX,), 


(4.22) 


+ ., 1 ™ |2 


We estimate the first term on the right hand side of the inequality inserting (4.21), 


lA+Xr < 2 


1 


ii/'" n+Ti'" -1- r)“T7"* -1- D+n”’ 


1 


1 


+ m 
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+ -1/2^2- [Wni/2 

+ 2|n™A^H^™+n-G(p™)|' 

< 4 


77-Dr 
2 2 


K, + l/2l^2+<, 


li/™ n+Ti’" -L li/™ n~n'^ j- ^ r)“ 


-'■1+1/2 




+ 4 




2'^ij + l/2-D2 + 2'*^i,i-l/2-D2 + 2 "^2 ki,i + l/2 1-^2 j 

+ 2|R™A,W^™+n™G(p™)|" 

< 8h^i/2,,-Di+n-IV 8|«ri/2,,^r</+ 

+ 8|<,_i/2-D2-n-/ + 2|n-A,M^- + n^,G{p^^)\^ 

< 8niax|V,W/-|{|R-i/2j + 111^1/2,, I l-Dr</ 

‘^i3 

+ K,+i/2l l-D2+n™f + 1-02”</} 

+ 2|R™A,T4^™+n™G(p™)|" 

Thus if we assume that At satisfies the CFL-condition (4.10), we have 

At^ lA+n-f < hY,{\uT+,n.\\Dt<,? + l<,+i/2ll^2+<,r} 

ij 

+ /i^|n-A,lT-+n-G(p-)|' 

Now summing (4.22) over all i,j, multiplying with and using the latter inequality, 
we obtain 

h^D+Y./::, = (l^-i/2oll^r<7l' + Ko-iWD^nZ?) 

+ h'^MY \DtnZ? + + 2G(p™ )) 

ij' 

+/i^+i E +<iG(pri)i' 

<C, 




where C > 0 is a constant independent of /i, thanks to the L^-bounds on and 
^hWh obtained in Lemma 4.1 and 4.3. This implies that 

Nt 

/i^'+^At EE(i '^T-1/2, + K,-i/2ll^2-<,P) < 

771=0 

Nt 

h'^M^YY.\DtN:,? <c. 

m—O i,j 

and therefore using Holder’s inequality and the uniform L°°-bounds on rih, (4.20). 
Using summation by parts, we realize that the other terms, (4.11) - (4.15) are in 
L°°([0,T];W~^'‘^{fl)) for q G [1,2*) where 2* = 2d/{d — 2) if d > 3 and any finite 
number greater than one if d = 2. □ 
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Remark 4.6. The preceeding lemma implies that the forward time difference of the 
approximation of the pressure D^ph = D^\nh\'^ is of the form D^ph = gn + kn 
where gh G L^([0,T] x O) and kh € L°°([0, T]; for any 1 < g < oo if 

d = 2 and for 1 < g < g* = 2d/{d — 2) if d > 2, uniformly in h > 0. Using this, we 
have that Wh = Uh + Vh where Uh and 14 solve 

-p.AhUh + Uh= gh, and - pAhVh + Vt = kh- 

By Lemma B.l, we have Uh,'^hUh G L^([0,T];L'^(U)) for 1 < q < d/{d— 1) 
and by standard results, 14, V?il4 G L°°([0, T]; L^(r2)). Hence DtWh,DtVhWh & 
Li([0,r];L9(H)) + L-([0,r];L2(r!)). 

Remark 4.7 (CFL-condition). The estimates from Lemma 4.3 imply that the ve¬ 
locity Uh := ^hWh G L°°([0, r]; (H)) uniformly in h > 0, 2* = 2d/{d— 2) or 

any number in [1, oo) if d = 2, using the Sobolev embedding theorem. Using an 
inverse inequality, we can bound it in the L°°{{0,T) x U)-norm as follows: 

max \uh\ < Ch~^ i / da: | < Ch~^ 

(rr,t)G(0.T)xn' \Jn J 

Thus the time step size At is of order ). In practice a linear CFL-condition 

seems to work well though. 


4.4. Passing to the limit d —^ 0. The estimates of the previous (sub)sections 
allow us to pass to the limit d —0 in a subsequence still denoted d, 

Uh ^ n > 0, inL'^([0,T] x fl), 1 < q < oo, 

Ph^P>0, inL'^([0, T] X U), 1 < q < oo, 

where ph ■= and 0 < n,p G L°°{[0,T] x fl). Using the “discretized” Aubin-Lions 
lemma A.l for 144 and 'S/hWh, we obtain strong convergence of a subsequence in 
L'J([0, T] X U) for any q G [0, oo) in the case of 144 and 1 < g < 2* in the case of 
hWh (2* = 2d/(d— 2) if d > 3 and any finite number greater than or equal to one 
if d = 2), to limit functions 14^, V14^ G L'^([0, T] x U). Moreover, from the estimates 
in Lemma 4.3 we obtain that 14^ G L°°([0,T] x U) n L°°([0,T]; Hence we 

have that (n, 14^,p) satisfy for any G C'^([0,T] x H), 

n Tvpt — n\7W ■ \7ipdxdt = — / / nG{p)<fdxdt 

i Jo Jn 

n Wip + pbS/W-Wii) dxdt = f ( pijjdxdt 

t Jo Jn 


where nG{p) is the weak limit of nhG{ph)- To conclude that the limit (n, 14^,p) is 
a weak solution of (1-6), we proceed as in the previous Section 4 and show that Uh 
in fact converges strongly: First, we recall that the limit n satisfies (3.9). 

On the other hand, from (4.22), we obtain (under the CFL-condtion (4.10)) 


DtK,^ < {|n”f + 

+ \d* {|< 




m |2 

J'l 


bj-il 




d2 

- 

2 2 


m I rA+_m |2 


2 

d2 

+ ^Dt + 4 ^ 2 + - 1 / 21 ^ 2 -<.] 


(4.23) 


+ |n-|^A,14^™+2|n-rG(p-) 
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Considering this inequality in terms of the piecewise constant functions Uh, Wh and 
Ph, multiplying it with a nonnegative C^-test function ip, integrating and then pass¬ 
ing to the limit /i —>■ 0, we obtain (using the bounds (4.20), the weak convergence 
of Uh and ph and the strong convergence of Wh and V^kC/j), 


n?ipt — • Vpdxdt < 


0 JQ, 


v?l\W -\- 2'n?G{p) I (fidxdt, (4.24) 


where denotes the weak limit of and n'^AW and •n?G{p) are the weak limits 
of n\AhWh and n\G{ph) respectively. 

Adding (3.9) and (4.24), we have 

— f f (ji'^ — n^'j (fit — {n^ — \/W ■ Vpdxdt 

< J J [2n?G{p) — 2nnG{p) + n'^AW — n^AW^ pdxdt. 

We now choose smooth test functions approximating p(t,x) = l[o,T](i)) where 
T G (0, T], in this inequality and then pass to the limit e —>■ 0 to obtain 

J {T)dx — J (^n'^{0,x) — n^{0,x)^dx 

<JJ (2n^G{p) — 2nnG{p) + AW + n'^AW — n^AW^dxdt. 

(4.25) 

By convexity of f{x) = x'^, we have > n^, on the other hand, the discrete 
L^-entropy inequality, (4.23), implies 

[ \nh{T,x)\^ dx < f \nl\'^dx+ f f {\nh\'^ AhWt + 2\nhfG{ph)) dxdt, 

t/O J0 J^ 

which gives, passing to the limit h —>■ 0, 

[ \n\‘^(T,x) dx < [ \no\^dx+ [ f f |npAW -I- 2|np( 


^Gip)) 


dxdt. 


Letting r —> 0, the second term on the right hand side vanishes (as the integrand 
is bounded), and we obtain 

f \n\'^{0,x) dx < ( \nQ\^dx 
JQ Jfl 

We deduce that |nP(0,-) = |rioP almost everywhere and that therefore the second 
term on the left hand side of (4.25) is zero. We have already estimated the first 
two terms on the right hand side of (4.25) in (3.13) and (3.14). To bound the other 
term, we use a discretized version of Lemma 3.7: 

Lemma 4.8. The weak limits {n,W,p) of the sequences {{'nh,Wh,Phy\h>o satisfy 
for any smooth function S' : K —>■ M, 

J (s(n)AW - dx=^J dx (4.26) 


where S{n)AW, S{n), pS{n) are the weak limits of S{nh)AhWh, S{nh) andphS{nh) 
respectively. 

Applying this lemma to the last term in (3.12) with S{n) = we can estimate 
it by 


/o JQ. 


n^AW — n^AW ] dx = — 
d 


— [ (pn'^ — pn'^'] dxdt 
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1 


< 0 , 




dxdt 


using again that by Exercise 3.37 in [24], rH v? < . Thus, 

J — n^^{T)dx < ^ 2 q ; + J J {n^ — dx dt. 

Gronwall’s inequality thus implies 

iTp- — {T)dx < 0 

By convexity of the function f{x) = x^ we also have almost everywhere 

and hence 

almost everywhere in (0, T) x O. Therefore we conclude that the functions Uh 
converge strongly to n almost everywhere, thus also p = v? and so the limit (n, W, p) 
is a weak solution of the equations (1.6). 



Proof of Lemma 4.8. We multiply the equation for Wh by S(nfi) and integrate it 
over the spatial domain 12, 


/ pAhWhS{nh)-WhS{nh)dx = - phS{nh)dx. 

Jn Ja 

Passing to the limit /i —> 0 in the last equation, we obtain 

I p.AWS{n) — WS{n)dx = — f pS{n)dx. 


(4.27) 


JQ JQ 

On the other hand, using [S{nfi) where ifs is a smooth mollifier converging 

to a Dirac measure at zero when 5 is sent to zero, as a test function in the weak 
formulation of the limit equation 


-pAW + W =p, 

and passing first to the limit d —>■ 0 and then h —>■ 0, we obtain 


/ pAWS{n)— WS{n) dx = — / pS{n)dx 
JQ JQ 

Combining the last identity with (4.27), we obtain (4.26). 


□ 


5. Numerical examples 

To test the scheme in practice, we compute approximations for the following two 
examples. 

5.1. Gaussian initial data. As a first example, we consider the initial data 

no{x) = ^(ixp{-10{xl+xl)) , (5.1) 

on the domain D = [—2.5,2.5]^ and h = 1/64 with pressure law p = v? and 
G{p) = 1 — p and p = 1. Strictly speaking, these are not homogeneous Neumann 
boundary conditions, but since the gradient of no near the boundary is very small, 
this works well in practice. 

In Figure 1 we show the approximations at times t = 0,1,2,4. We observe 
that the cell density in the middle hrst reaches the maximum possible and then 
starts spreading with a relatively narrow transition region between zero density 
and maximum density. 
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t = 0 t = 1 



Figure 1 . The approximations of the cell density n for initial data 
(5.1) on n = [—2.5, 2.5]^ with mesh width h = 1/64. 


5.2. Two Gaussians. As a second example, we use the inital data consisting of 
two Gaussian pulses with centers at a; = (0.7,0) and x = (—0.6,0.2), 


no{x) 


i exp (-10 ((xi - 0.7)^ + X 2 )) 

+ i exp (-20 ((xi + 0 . 6 )^ + (x 2 - 0 . 2 )^)) 


(5.2) 


on the same domain, fl = [—2.5,2.5]^, with /r = 1, pressure law p = and 
G{p) = 1—p and mesh width h = 1/64. The approximations computed at times 
t = 0, 2,4 ,6 are shown in Figure 2. The interface between the area with maximum 
cell density and zero cell density seems to be sharper than in the previous example, 
this appears to be caused by the pressure law with the higher exponent 7 . Further 
tests with higher and lower exponents confirmed that assertion. 


Appendix A. Discretized Aubin-Lions lemma 

Lemma A.l. Let Uh '■ Lt be a piecewise constant function defined on a grid 

on [0, T) X LI, Q a bounded rectangular domain, satisfying 

pT r 


/ / \uh\'^ + \VhUh\'^dxdt<C 


(A.1) 


Jo Jn 

for some 00 > q > 1, uniformly with respect to h > 0 and 


DtUh — -Afifh + <?/i + kh, (-'^•2) 

where Ah is a first order linear finite difference operator, and fh,gh,kh ■ Lt ^ 
are piecewise constant functions, satisfying uniformly in h > 0, 
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Figure 2. The approximations of the cell density n for initial data 
(5.2) on n = [—2.5, 2.5]^ with mesh width /i = 1/64. 


[ [ + \9hP + \kh\dxdt < C, 

Jo Jci 

for some oo > ri, r 2 > 1. Then Uh ^ u in ^^*^([0, T) x il). 


(A.3) 


Proof. Denote uu a piecewise linear interpolation of Uh in space piecewise constant 
in time and similarly, let gh, fh and kh piecewise linear interpolations of gh, fh and 
kh respectively in space and piecewise constant in time such that 

DtUh = Ah fh +gh + kh. (A.4) 


By Ladyshenskaya’s norm equivalences [21, p. 230 ff], we have 


\uh\\wi.g(^Q)dt < C j J \uh\‘‘+ hUh\‘^ dxdt 


2’■l{Q) + \\9h\\L-2(n) + \\f^h\\L^n)dt<C / {hr + \ghp + \kh\dxdt 


0 Jn 


where the right hand sides are bounded by assumptions (A.l) and (A.3). Since 
L^{n) C 1F“^’®(D) for 1 < s < 1* = d/{d—l), we have that kh € T^([0, T]; 1F“^’®(D)) 
for 1 < s < 1* = (i/(d — 1) and hence thanks to this and (A.4), we obtain 

Uh G L«([0,T);lFi-«(D)), DtUh G V([0, T); (H)), 


uniformly with respect to the discretization parameter h > 0. Thus we can ap¬ 
ply the version [11, Theorem 1] of the Aubin-Lions lemma to find that up to a 
subsequence —>■ m in L‘^([0,T) x 11) and the limit u G L'^{[0,T);W^'‘>{i})). By 
[21, Lemma 3.2., p. 226] this implies that also Uh ^ u in L'?([0,T) x 11) (and 
VhUh ^ Vu ). □ 
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Remark A.2 (Derivatives). If the Uh in Lemma A.l are of the form VhVh for some 
Vh piecewise constant function, this lemma implies that VhVh Vu in L‘>, again 
applying [21, Lemma 3.2., p. 226] 

Appendix B. Technical Lemmas 


In this section, we prove the following lemma: 

Lemma B.l. Let Uh solve the difference equation 

- (LwhiAhVhUh) + ChUh = fh, X G D, (B.l) 

with homogeneous Neumann boundary conditions, where Ah is a diagonal positive 
definite d x d-matrix with entries > rj > 0 and c^, > u > 0 uniformly in h > 0, 
X G LI is a rectangular domain in and 

WfhWL^n) < M, 

uniformly in h > 0. We have denoted V/j := Nf and div^i := div)[( (or alternatively 
Vh := and div^ := divf). Then 

||wii||L‘i(n) + IIVhM/t|U‘!(n) < C, 

where 1 < q < d/{d — 1), for a constant C > 0 independent of h > 0. 

The proof of this lemma will be a (simplified) finite difference version of the 
proof of Theorem 2.1 in [4]. But before proving the lemma, we need to introduce 
some notation. 


Notation B.2. For any r G (l,oo), we denote by the Marcinkiewicz space 

with norm defined by 

llM||L-.~(n) = supA|{x e D : \u{x) > A}|^/’'. 

A>0 

The Marcinkiewicz spaces are continuously embedded in L^{Ll) for any 1 < 9 < r, 
[15]: 

l|M||L9(n) < C'(g,r, |D|)||M||ir.,=o(Q), qG[l,r). (B.2) 

Moreover, we need the trunctation operator Sk defined as follows: 


Notation B.3. Let fc > 0 be a real number. Then we define the truncation operator 
5*; : K —>■ K by 


Sk{s) = 


|s| ’ 


if |s| < k, 
if jsj > k. 


It will be convenient in the proof to use the following tuple notation for the finite 
difference approximations: 


Notation B.4. We denote i := (ii,..., id), it = 1, ■.., Ng, Ni the number of cells 
in the .^th spatial direction, a d-dimensional tuple and and Ui the approximation in 
cell Ci := ((fi — \)h,iih] x • • • x [{id — \)h,idh]. The piecewise constant function 
Uh can be written as 

Uh{x) :='^Uilcfix), xGLL. 

i 

We also need the following auxilary result: 


Lemma B.5. Let Uh solve the difference equation (B.l) under the assumptions of 
Lemma B.l. Then 


f \VhSk{uh)f+ \Sk{uh)fdx<CMk, Vfc>0, 
Jn 

for some constant C > 0 independent ofh>0. 


(B.3) 
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Proof. Given fc > 0, we multiply equation (B.l) by Sk{uh) and integrate over the 
domain Q. After changing variables in the integrals, we obtain 

/ {AhVhUh)-'^hSk{'u.h) + ChUhSk{uh)dx= / fhSk{uh)dx. (B.4) 
Jo. JQ. 

The right hand side can be bounded hy Mk using Holder’s inequality. The left 
hand side, we can rewrite and estimate as follows 

/ hUh) ■ hSk{uh) + ChUhSk{uh) dx 

Jn 

= / {AhVhSk{uh))-yhSkiuh)+ Ch\Sk{uh)\'^ dx 
JQ 

+ / {Ah{\^h[uh-Sk{uh)]))-'^hSk{uh) + Ch{uh-Sk{uh))Sk{uh)dx 
JQ 

> hSk{uh)\\ 'l'^(Q) + ^\\Sk{uh)\\'L2(Q) 

+ / {AhC^h[uh-Sk{uh)]))-'^hSk{uh) + Ch{uh-Sk{uh))Sk{uh)dx. 

JQ 

(uh — Sk{uh)) is either zero or has the same sign as Sk{uh). Therefore (ut — 
Skiuh))Skiuh) > 0 and 

/ Ch{uh-Sk{uh))Sk{uh)dx>0. 

JQ 

In order to prove that the other term is positive as well, we will show that 

DJSk{ui)DJ {u^-Sk{u^)) > 0, Vi, i=l,...,d. 

The proof of this fact consists of boring case distinctions and is exactly analoguous 
for ^ = 1, 2, (3), therefore we will do it only for £ = 1 and omit writing the tuple 
index i. Then we have 


{ui - Skiui))D^ Skiui) = < 


{u. 

i - k){k - Ui-i), 

Ui 

> k, |uj_ 

i| < k, 

{u, 

i + k){-k - 

Ui 

-.if 

1 

V 

-il < k, 

0, 



;| < fc, \u^. 

-i| < k, 

(- 

Ui-i + k){ui - k), 

\u^ 

; < k, Ui- 

-1 > k, 

{- 

Ui-i - k){ui + k), 


\ < k, Ui- 

1 

V 

0, 


Ui 

> k, Ui-i 

> k, 

0, 


Ui 

1 

V 

1 

V 

(M: 

i - Ui-i - 2k)2k, 

Ui 

> k, Ui-i 

1 

V 

-( 

Ui — Ui-i + 2k)2k, 

Ui 

< —k, Ui- 

-1 > k. 


The potential reader is welcome to check that these are all the possible cases and 
that each of the terms on the right hand side is nonnegative. Thus we have that 

/ {AhVhUh)-^hSk{uh) + ChUhSk{uh)dx 
JQ 

> r]\\'^hSk{uh)\\‘j^2(^Q'i + k'\\Sk{uh)\\L2(^fi^ 

which implies (B.3) together with the estimate on the right hand side of (B.4) □ 

Proof of Lemma B.l. First, we note that by the discrete Gagliardo-Nirenberg-Sobolev 
inequality, [2, Thm. 3.4], 


\Skiuh)\‘^ dx < 


\'^ hSk{uh)\'^ + \Skiuh)\'^dx 
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where 2* = 2d/{d — 2) if d > 3 and any number with 1 < 2* < oo if d = 2, and 
where C is a constant depending on |r2| but not on d > 0. By Lemma B.5, we can 
bound the right hand side and obtain therefore 

/ |5feK)P*dx<C(fcM)^. (B.5) 

Jn 

Now we define the set B{k) by 

B{k) = {CiCn: |wi| > k}. 

We have 

[ \Skiuh)f dx > k^'\B{k)\, 

JBik) 

and therefore, using (B.5), 

\B{k)\ < / \Sk{uh)\‘^ dx < [ \Skiuh)\‘^ dx < (B.6) 

JBik) ^ Ja k 2 

which implies that Uh G for r = 2*/2 (which is d/(d — 2) if d > 3) since 

the choice of A: > 0 was arbitrary. Now denote 

as(fc) := {q C : 3 j, li - JI = 1, \u^\ > k} 

WH) :=B{k)UdB{k), 

B{kY := n\B{k), 

where \i — j\ = maxi<£<£i \ii — ji\. Informally speaking, the cells in dB{k) have a 
neighbor cell which is contained in B{k). We have 

CMt 

\dB{k)\<{i'^-l)\B{k)\<—^, 

k 2 

by (B.6). Now let A > 0, fc > 0 and decompose 


{a; € 12 : \\IhUh{x)\ > A} = {a; G 12 : \\IhUh{x)\ > Aanda: G B{k)} 

U {a; G n : \VhUh{x)\ > Aandx G B{kY}. 

Hence 


|{a; G H : |V;iUft,(a;)| > A}| < \B{k)\ + |{x G 12 : \\7 hUh{x)\ > Aanda; G B{kY}\. 

On B{kY and the cells bordering the set, we have \uh\ < k and therefore Uh = 
|<5'fe(u/i)|- Hence we can estimate the size of the second set in the above inequality. 


|{a; G 12 : hUh{x)\ > Aanda; G B{kY}\ 

= |{a; G S2 : \'^hSk{uh){x)\ > Aandx G B{kY}\ 

< |{a; G 12 : \VhSk{uh){x)\ > A}| 

- ^ J \^hSk{uh)\'^dx, 

where we have used Chebyshev inequality for the last step. Now we can estimate 
the size of the set {a; G 12 : \\7hUh{x)\ > A} using (B.3) once more, 


|{a; G 12 : |V/jWft(a;)| > A}| < 


k^ 


Choosing k = A^^+s, we obtain 


CkM 


A^fS |{a; g : |V/tM/i(a:)| > A}| < C{d,M, |12|). 


If d > 3, we have and so Uh, VkUt G L’'’°°(12) for 1 < r < d/(d — 1). 

For d = 2, since 2* is an arbitrary finite positive number, we can achieve the same. 
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Using the embedding of the Marcinkiewicz spaces, (B.2), we obtain the claim of 
the lemma. □ 
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